Spatial optical solitons in nonlinear photonic crystals 
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We study spatial optical solitons in a one-dimensional nonlinear photonic crystal created by an 
array of thin-film nonlinear waveguides, the so-called Dirac-comb nonlinear lattice. We analyze 
modulational instability of the extended Bloch-wave modes and also investigate the existence and 
stability of bright, dark, and "twisted" spatially localized modes in such periodic structures. Addi- 
tionally, we discuss both similarities and differences of our general results with the simplified models 
of nonlinear periodic media described by the discrete nonlinear Schrodinger equation, derived in the 
tight-binding approximation, and the coupled-mode theory, valid for shallow periodic modulations 
of the optical refractive index. 



I. INTRODUCTION 



Discrete spatial optical solitons have been introduced 
and studied theoretically as spatiaUy locahzed modes of 
periodic optical structures (see, e.g., Refs. and also 
a review paper |^), and they have recently been ob- 
served experimentally in arrays of nonlinear single-mode 
optical waveguides |6|. A standard theoretical approach 
in the study of the discrete spatial optical solitons is 
based on the derivation of an effective discrete nonlin- 
ear Schrodinger (DNLS) equation and the analysis 
of its stationary localized solutions - discrete localized 
modes In the solid-state physics, the similar approach 
is known as the tight-binding approximation which, in ap- 
plication to optical waveguide arrays, corresponds to the 
case of weakly coupled fundamental modes excited in each 
waveguide of the array. The analogous concepts appear 
in other fields such as the nonlinear dynamics of the Bose- 
Einstein condensates in optical lattices ||^. 

On the other hand, weak nonlinear effects in optical 
fibers with a periodic modulation of the refractive index 
(often called optical grating) are well studied in the frame- 
work of another approach, the coupled-mode theory. The 
coupled-mode theory is based on a decomposition of the 
electric field into the forward and backward propagating 
components, under the condition of the Bragg resonance. 
Such an approach is usually applied to analyze nonlinear 
localized waves in the systems with a weakly modulated 
optical refractive index known as gap (or Bragg) soli- 
tons 1^, and such gap solitons are known to appear in 
other fields §. 

Thus, the theory of spatial and temporal optical soli- 
tons in periodic structures developed so far is based on 
one of the two approaches, the DNLS equation or the 
coupled-mode theory. However, real experiments in the 
nonlinear guided-wave optics are conducted in the peri- 
odic structures of more complicated geometries and un- 
der the conditions when none of those approximations are 
valid. In such a case the applicability of the tight-binding 
approach and the corresponding discrete equations, from 



one hand, and the coupled-mode theory, from the other 
hand, become questionable, especially for the analysis of 
the linear stability of nonlinear localized modes. There- 
fore, a consistent theory of nonlinear effects and localized 
modes in periodic media is still missing. 

One of the main features of wave propagation in pe- 
riodic structures (which follows from the Floquet-Bloch 
theory) is the existence of a set of forbidden band gaps in 
the transmission spectrum. Therefore, the nonlinearly- 
induced wave localization can become possible in each of 
these gaps. However, the effective DNLS equation de- 
rived in the tight-binding approximation describes only 
one transmission band surrounded by two semi-infinite 
band gaps and, therefore, a fine structure of the band- 
gap spectrum associated with the wave transmission in a 
periodic medium is lost. On the other hand, the coupled- 
mode theory of gap solitons |^ describes only the modes 
localized in an isolated narrow gap, and it does not allow 
to consider simultaneously the gap modes and conven- 
tional guided waves localized due to the total internal re- 
fiection. The complete band-gap structure of the trans- 
mission spectrum and simultaneous existence of local- 
ized modes of different types are very important issues in 
the analysis of stability of nonlinear localized modes [0 . 
Such an analysis is especially important for the theory of 
nonlinear localized modes and nonlinear waveguides in 
realistic models of nonhnear photonic crystals (see, e.g., 
the recent paper and references therein). 

In this paper, we consider a simple model of nonlinear 
periodic layered media where a periodic optical structure 
is formed by an array of thin-film nonlinear waveguides 
embedded into an otherwise linear dielectric medium (see 
also Ref. [0). Such a structure can be regarded as a 
nonlinear analog of the so-called Dirac comb lattice , 
where the effects of the linear periodicity and band-gap 
spectrum are taken into account explicitly, whereas non- 
linear effects enter the corresponding matching condi- 
tions allowing a direct analytical study. 

We analyze nonlinear localized modes in an infinite 
structure consisting of a periodic array of nonlinear 
waveguides, similar to the geometry of the experiments 
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with discrete optical solitons Q. First, we study modula- 
tional instability of extended modes in both self-focusing 
and self-defocusing regimes. Then, we discuss different 
types of nonlinear localized modes (such as bright, dark 
and "twisted" spatial solitons) and also analyze numer- 
ically their linear stability. We emphasize both similar- 
ities and differences between our results and the results 
obtained in the framework of the DNLS equation and the 
continuous coupled-mode theory. 



II. GENERAL APPROACH 



A. Model 



and field amplitude, respectively. Then, the normalized 
nonlinear equation has the form 



.dip 92^/; , - 



(2) 



where the real function T{I\x) = (P'D~^[£{X) — £ -I- 
g{X)I\EQ\'^] describes both nonlinear and periodic prop- 
erties of the layered medium, and / = is the normal- 
ized local wave intensity. We note that the system (^) is 
Hamiltonian, and for spatially localized solutions it con- 
serves the total power. 



P = 



-t-CXD 



|'(/'(a;, z)p dx. 



We consider the electromagnetic waves propagating 
along the Z-direction of a slab-waveguide structure cre- 
ated by a periodic array of thin-film nonlinear waveg- 
uides (see Fig. |^). Assuming that the field structure in 
the Y direction is defined by the linear guided mode of 
the slab waveguide £{Y; X), we separate the dimensions 
presenting the electric field as E{X,Z)£{Y\X). Then, 
the evolution of the complex field envelope E{X^ Z) is 
governed by the nonlinear Schrodinger (NLS) equation. 



.dE d^E 



eiX)E + g{X)\E\^E = 0, 



(1) 



where D is the diffraction coefficient {D > 0) . The phase 
velocity of the guided waves is defined by the function 
£{X), whereas g{X) characterizes the Kerr- type nonlin- 
ear response of the layers. We assume that either the 
function e{X) or g{X) (or both of them) is periodic in 
X, i.e. it describes the periodic layered structure simi- 
lar to the so-called transverse Bragg waveguides created 
by the nonlinear thin-film multilayer structures |14|,|l5 | or 
the impurity band in a deep photonic band gap 




FIG. 1. Array of thin-film nonlinear waveguides embed- 
ded into a linear slab waveguide. Gray shading shows a profile 
of a nonfinear localized mode. 

In order to reduce the number of physical param- 
eters, we normalize Eqs. (|^) as follows: E{X, Z) = 
ipix, z)Eoe^ '^^ , where £ is the mean value of the function 
s{X), X = X/d and z = ZD/d^ are the dimensionless co- 
ordinates, d and Eq are the characteristic transverse scale 



At this point, it is important to mention that Eq. (Q) 
describes the beam evolution in the framework of the so- 
called parabolic approximation, valid for the waves prop- 
agating mainly along the z direction (see also Ref. [TtI ] 
and discussions therein). In other words, the characteris- 
tic length of the beam distortion due to both diffraction 
and refraction along the z axis should be much larger 
than the beam width in the transverse direction x. This 
leads to the condition of a weakly modulated periodicity, 
|e(X)-£|«|£|. 

We look for stationary localized solutions of the nor- 
malized equation (0) in the standard form 



ip{x,z) =u{x;P)e'^\ 



(3) 



where p is the propagation constant, and the amplitude 
function u(x] (3) satisfies the stationary nonlinear equa- 
tion: 



d^u 

(3u + + T{I]x)u = 0. 
dx'^ 



(4) 



If there is no energy flow along the transverse direction 
X, then the function u{x) is real, up to a constant phase 
which can be removed by a coordinate shift z ^ z — zq. 
This is always the case for spatially localized solutions 
with vanishing asymptotics, u{x — > ±00) — 0. 

To simplify our analysis further, we assume that the 
linear periodicity is associated only with the presence of 
an array of the thin-film waveguides, and define the re- 
sponse function in the model Eq. (0) as follows 



T{I;x)^ {a + -fI)Six- hn), 



(5) 



where h is the spacing between the neighboring thin-film 
waveguides (the lattice period), and n is integer. The 
total response of the thin-film layers is approximated 
by the delta-functions, and the real parameters a and 
7 describe both linear and nonlinear properties of the 
layer, respectively. Without loss of generality, the non- 
linear coefficient 7 can be normalized to unity, so that 



2 



7 = +1 corresponds to self-focusing and 7 = — 1 to self- 
defocusing nonlinearity. The linear coefficient (a > 0) de- 
fines the low-intensity response, and it characterizes the 
corresponding coupling strength between the waveguides. 
Model (^),(^ can be regarded as a nonlinear analog of 
the Dirac comb lattice, earlier studied in application to 
the photonic crystals in the linear regime only pSl . 



B. Dispersion Properties and Discrete Equations 

Following the path outlined in Ref. |f 7]| , we present 
the stationary modes defined by Eqs. (j^M^ as a sum of 
the counter-propagating waves in each of the linear slab 
waveguides. 



u{x) = a„e-''(^-"'') + 6„e+^(^-"''), 



(6) 



where nh < x < (n + l)h. Then, we express the coeffi- 
cients a„ and 5„ in terms of the wave amplitudes at the 
nonlinear layers, u„ = u{hn), 



2 sinh (fih) 

^71 ; 



(7) 



where fi{(3) = Finally, we substitute Eqs. (||),(0) 

into Eqs. and find that the normalized ampli- 

tudes Un = \/\^l\un satisfy a stationary form of the 
DNLS equation: 

TjUn + [Un-l + Un+l) + xWn^U,, = 0, (8) 

where x — sign(^7), and 

■q{(3) = -2 cosh(^/i) + aC(/9), 

(9) 

^{(3) = sinh(/i/i)//i. 




FIG. 2. Characteristic dependencies of the parameters 77 
and f on the propagation constant /3. Gray shading marks 
the transmission bands. The lattice parameters are h = 0.5 
and a = 10. 



Linear solutions of Eq. (when the term ~ xl^^nP 
vanishes) have the form C/„ = C/oe'^", where K = 
±cos~^(— 7//2) is the wave number. Therefore, extended 
linear solutions with real K can exist for \rj\ < 2. On 
the other hand, nonlinear localized modes with exponen- 
tially decaying asymptotics can appear only for I77I > 2 
(when K is imaginary), this condition defines the band- 
gap structure of the spectrum. Characteristic dependen- 
cies of rj and ^ vs. the propagation constant /3 are pre- 
sented in Fig. ^, where the bands are shown by a gray 
shading. The first (semi-infinite) band gap corresponds 
to the total internal reflection (IR). On the other hand, 
at smaller /3 the spectrum band gaps appear due to the 
resonant Bragg-type reflection (BR) from the periodic 
structure. 



C. Linear Stability Analysis 

To study linear stability of localized modes, we should 
consider the evolution of small-amplitude perturbations 
of the localized state presenting the solution in the form 



ij{x,t) = {^u{x) + v{x)e'^' +w*ix)e~'^''Y 



if3z 



(10) 



From Eqs. (^ and (^, we obtain the linear eigenvalue 
problem for small v{x) and w{x): 

+ 00 



[(a + 27|u„|^)u + 7u^w] (5(x - /in) = 0, 



(11) 



+ 00 



^ [(a -F 2-f\un\^)w -f 7«)'v] 6ix - hn) = 0. 



After representing the fields as sums of the counter- 
propagating waves (see Sec. II B), we reduce Eqs. ( [Tl| ) 
to a set of the discrete equations for the amplitudes at 
the layers: 



77(/3 -I- T)Vn + {Vn-l + Vn+l) + 

+7C(/3 + r) [2|u„|2w„ + ulwn] = 0, 

T]{I3 - T)Wn + [Wn-l + W„+l) + 

+7C(/3 - r) [2|u„|2u;„ + [KYvn] = 0. 



(12) 



In general, the solutions of this eigenvalue problem fall 
into one of the following categories: (i) internal modes 
with real eigenvalues that describe periodic oscillations 
("breathing") of the localized state, (ii) instability modes 
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that correspond to purely imaginary eigenvalues with 
ImF < 0, and (iii) oscillatory unstable modes that appear 
when the eigenvalues are complex (and ImF < 0). Addi- 
tionally, there can exist decaying modes when ImF > 0. 
However, from the structure of Eqs. (|l^) it follows that 
the eigenvalue spectrum is invariant with respect to the 
transformation F — ±F*, if all the amplitudes Mne**^ are 
real, where 4> is an arbitrary constant phase. In such 
a case exponentially growing and decaying modes always 
coexist, and the latter do not significantly affect the wave 
dynamics. 

It is important to note that the grating properties 
are expressed through the functions r]{i3) and ^(/3) [see 
Eq. dgj)] in both stationary [Eq. (||)] and perturbation 
[Eq. equations. Therefore, the functions r]{P) and 

^(/3) fully characterize the existence and linear stability 
properties of localized and extended solutions of the orig- 
inal model (||),(||) with stationary intensity profiles. 



III. APPROXIMATE MODELS 
A. Tight-Binding Approximation 

When each of the thin-film waveguides of the periodic 
structure supports a fundamental mode that weakly over- 
laps with the similar mode of the neighboring waveguide, 
the modes become weakly coupled via a small change of 
the refractive index in the waveguides. Then, the mode 
properties can be analyzed in the framework of the so- 
called tight-binding approximation^ well developed for the 
problems of solid-state physics . 

In order to employ this approximation in our case, first 
we analyze the properties of a single thin-film waveguide 
in the linear regime and solve Eq. (^ with J- {I; x) = 
aS{x) to find the spatial profile of a linear guided mode, 

Us{x) — exp (— a|a;|/2) . 

and the corresponding value of its propagation constant, 
Ps = o? 1^. Second, we consider the interaction between 
the waveguides in the array assuming that the total field 
can be presented as a superposition of slightly perturbed 
waves localized at the isolated waveguides. Specifically, 
we assume that the propagation constant remains close 
to its unperturbed value (i.e. |«c?V'/di + /3s'0l ^ I/^sV'Dj 
and neglect small variations in the spatial profiles of the 
localized modes. Then, we seek a general solution of 
Eqs. (§),(||) in the form 



ip{x,z)^ ^ ipn{z)us{x - nh). 



(13) 



n— — oo 



In such an approximation, the wave evolution is charac- 
terized by the amplitude functions ipniz) only. In order 
to find the corresponding evolutionary equations, we sub- 
stitute Eq. ( p^ ) into Eqs. (||) and (||) and, multiplying the 



resulting equation by Us{x — mh), integrate it over the 
transverse profile. According to the original assumption 
of weakly interacting waveguides (valid for ah 3> 1 and 
l7l|V'n(-z)P <C a), in the lowest-order approximation we 
put \ip{nh, z)\'^ ~ |'(/'n(-z)P and neglect the overlap inte- 
grals, Us{x — nh)us{x — mh)dx, with |n — m| > 1, 
which produce the terms of higher orders. Finally, we de- 
rive a system of coupled discrete equations for the field 
amplitudes at the nonlinear layers (up to small pertur- 
bations, since ^l){nh, z) — ipniz) + 0(e~"''^^)): 



+ Ps^n + ^e-"''/2(V.„_i + ^„+i) 
dz 2 



(14) 



Stationary solutions of Eq. (|lj) have the form ■;/;„ — 
UnC^^'^ , where the amplitudes u„ are given by Eq. (||) 
with 

'/ = -^^/3-^%-"''^^ e=-e-"'^/^ (15) 



Relations ( [l5|) can be found as a series expansion of the 
original dispersion relations (^) near the edges of the first 
transmission band, in the limit ah 3> 1. 



B. Coupled-Mode Theory 

Now we consider the opposite limit ah <^ 1, when the 
first Bragg-reflection gap is narrow, i.e. (/3i — 02) ^ 
\(3i,2\, where /3i^2 are the propagation constants at the 
gap edges, defined by the condition 77(/3i,2) = 2. From 



Eq. 



it follows that (32 



-{n/hy and pi P2 + 



2a/h + 0{a^^^). To find solutions close to the BR gap, 
we present the total field in the form 

ipix, z) = ai{x, z)ub{x; Pi) + 02(2;, z)ub{x; P2), (16) 

where Oj are unknown nonlinear amplitudes, and 
Ub{x]pi^2) are the linear Bloch functions which satisfy 
Eqs. (||) and (||) at 7 = 0. The Bloch functions can 
be found in an explicit form: Ub{x] P2) = sin^xn/h), and 
Ub{x + nh;Pi) = (-1)" sin[(V2-x) for < x < h. 

Note that the field amplitudes at the layers are u„ ~ 
(— l)"oi(n/i), since Ub{nh; P2) = 0. 

In order to find the equations for the amplitudes 
aj{x,z], we substitute Eq. ( [l^ ) into the original model 
Eqs. (|^) and (^. Next, we use the fact that the gap is 
narrow, and close to its edges the Bloch functions are 
weakly modulated, i.e. \dai^2/dx\ <C \ai,2/h\. This as- 
sumption allows us to keep only the lowest-order terms. 
Then, we multiply the resulting equation by Ub{x;Pi^2) 
and integrate it over one grating period. Finally, the 
coupled equations for the modulation amplitudes are 



.dai 27r 9a2 2 ^ 

.802 , „ 2Trdai 

(- P2a2 — r^r~ = 0' 

az li ax 



0, 



(17) 
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Equations ( p7| ) derived above allow a direct compar- 
ison between the coupled-mode theory and the general 
results for the stationary localized solutions for which 
ai_2(a;,z) — a\j,{x)e^f^^ . Additionally, since the func- 
tions aj are weakly modulated, the spatial derivatives 
can be approximated by finite differences between the 
amplitudes at the layers (note that such an approxima- 
tion is only valid for /3 ~ /3i,2)- After simple algebra we 
obtain a discrete equation for the field amplitudes at the 
layers which has the form of the DNLS equation with 



77(/?) = 2-— (/3-/3i)(/3-/32), 



(18) 



Similar to the case of the tight-binding approximation, 
the dispersion relations (^|) can be found from the gen- 
eral result (^) by performing a series expansion near the 
band edge value /?2- 



C. Two-Component Discrete Model 

The principal limitations of both the tight-binding ap- 
proximation and the coupled-mode theory is explained by 
the fact that they are both valid in local narrow regions 
of the general band-gap structure and under special as- 
sumptions. Indeed, these two approaches are applicable 
when the dimensionless parameter ah is either small or 
large, and none of those approaches covers the interme- 
diate cases. While a general study should rely on the nu- 
merical solutions with the exact dispersion relations (S), 
it is useful to consider a simplified model which can (at 
least, qualitatively) describe the wave properties close to 
the BR gap (/3 — /3i) in the transitional region, being 
valid for ah ~ 1 as well. To achieve this goal, we ex- 
tend the tight-binding approximation and the correspond- 
ing DNLS equation introduced above in Sec. Ill A. 

We note that Eq. ( p^ ) can be considered as a rough dis- 
cretization of the original model (^) , with only one node 
per grating period located at x = nh. Then, the natural 
generalization is to include additional nodes located be- 
tween the nonlinear layers at the positions x = (n-\-l/2)h. 
Since the refractive indices at the node position are now 
different, we obtain a new system of coupled discrete 
equations which correspond to a two-component super- 
lattice: 



+ Mn + Pi(V'«-i/2 + V'n+1/2) 

+l\'ipn?i^n = 0, (19) 
+ /?2V'ri+l/2 + P2(V'« + ^n+l) = 0, 

We find that, similar to the general case, the stationary 
mode profiles can be expressed in terms of the amplitudes 



u„, which satisfy the normalized DNLS Eq. (|8|) with the 
following parameters: 



7/(/3) = 2-(pip2)-'(/?-/3i)(/3-/32), 



(20) 



One can immediately see that our new model ( [19|) de- 
scribes an effective system with a semi-infinite IR gap 
and a BR gap of a finite width. We note that, although 
dispersion relations in Eqs. (|2^) and ( p8| ) look similar, 
the latter relation is only valid for (3 ~ /3i,2, so that 
the coupled-mode theory describes only a single isolated 
Bragg-reflection gap. Therefore, the model (|l9| ) provides 
an important generalization of the DNLS theory, which 
also has a wider applicability than the coupled-mode the- 
ory. 

Because the model ( p^ describes an extra gap in the 
transmission spectrum, a comparison with the original 
model becomes more complicated. We choose the model 
parameters pj and 7 in order to match the dispersion re- 
lations in the vicinity of (3i , and satisfy a relation between 
Un and Un+i/2 following from Eqs. (^),(0): 



Pi 



-2 cosh{y^h / 2) /{dr]/dl3)p^ 



P2 = {f3i - /32)/[2cosh(VA/i/2)], 



(21) 



where the functions ry(/3) and ^(/3) should be calculated 
according to Eq. (^). 

The three different approximate models discussed 
above allow a simple analysis of the limiting cases, and 
they will be used below in calculating different properties 
of the extended and localized modes. 



IV. MODULATIONAL INSTABILITY 

A. General Analysis 

First, we analyze the properties of the simplest ex- 
tended (plane-wave) solutions of the model (||) and ^ 
which have equal intensities at the nonlinear layers, 
-^0 = = const, and correspond to the first trans- 

mission band. These solutions have the form of the 
so-called Bloch waves (BWs) Un = upe'^", where the 
wave number K is selected in the first Brillouin zone, 
\K\ < n. Using Eq. (||), we find the dispersion relation 
as if = ±cos^^[ri{(3;a)], where a = (a -I- 7/0) defines 
the layers response modified by nonlinearity. Since the 
transmis sion bands are defined by the condition |r/| < 2 
(see Sec. [IB ), the band structure shifts as intensity in- 
creases. Indeed, by resolving the dispersion relation we 
determine a relation between the propagation constant 
and the wave intensity. 



loW = - 



[2 cos K -^T]{f3)] 



(22) 
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Since in the first transmission band j3 > — (tt/Zi)^, we 
have > (see also Fig. ||), and it foUows from Eq. (||) 
that the propagation constant (3 increases at higher in- 
tensities in a self- focusing medium (7 > 0), and decreases 
in a self-defocusing medium (7 < 0). 

One of the main problems associated with the non- 
linear BW modes is their instability to periodic modu- 
lations of a certain wavelength, known as modulational 
instability (see also In order to describe the sta- 

bility properties of the periodic BW solutions, we anal- 
yse the evolution of weak perturbations described by the 
eigenvalue problem ([l^). Due to periodicity of the back- 
ground solution u{x), it follows from the Bloch theorem 
that the eigenmodes of Eq. ( [l^ ) should also be periodic, 
i.e. v{x + h)^ v{x)e''^i+^^ and w{x + h) = w{x)e''-'^-^^ . 
Then, we obtain the following solvability condition: 

[77(/3 + r) + 27^(/3 + r)/o + 2 cos(g + K)] 

X [r,iP - r) + 27e(/3 - r)/o + 2 cos(g - K)] (23) 

= 7'e(/3 + rK(/3-r)/2. 

Possible eigenvalues T are determined from the condi- 
tion that the spatial modulation frequencies q, which 
are found from Eq. (p3|), are all real. Therefore, the 
eigenvalue spectrum consists of bands, and the instability 
growth rate can only change continuously from zero to 
some maximum value. We also note that the spectrum 
possesses a symmetry F — > ±F*, and it is sufficient to 
study only the solutions with ReF > 0. 



Q(F) = cosb(F)], 
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FIG. 3. Modulation frequency q vs. real (left) and imagi- 
nary (right) parts of the perturbation eigenvalue F, for the 
staggered BW modes in a self- focusing medium (K = tt, 
/3 = 7, Q = 3, ft = 0.5, 7 = +!)■ Solid lines — stable 
(ImF = 0) and dashed — unstable (ImF 7^ 0). 



B. Stability of Staggered and Unstaggered Modes 

In what follows, we consider two characteristic cases 
of the stationary nonlinear BW modes, when they are 
(i) unstaggered {K = 0) or (ii) staggered (K — tt). To 
describe a transition from purely real to complex linear 
eigenvalues, we consider the function 



defined from Eq. (|23|), and extend the solution from the 
real axis to the complex plane by writing a series ex- 
pansion: g(ReF -I- ilmF) = Q(ReF) -h (5'(ReF)(iImF) + 
(l/2)Q"(ReF)(iImF)2 + 0(ImF3), where the prime de- 
notes differentiation with respect to the argument. Since 
the modulation frequency should remain real, the sec- 
ond term in the series expansion should vanish. Then, 
we conclude that complex eigenvalues only appear at the 
critical points, where 



Q'(F) - 0. 



(24) 



Therefore, (in)stability can be predicted by studying the 
function (3(F) on the real axis only, and then extending 
the solution to the complex plane at the critical points 
(if they are present) determined by Eq. (HJ). The real 
modulation frequencies are found as g = cos^^((3) in the 
interval — 1 < Q < 1, therefore MI only appears at the 
critical points where q'(T) = or g = 0, tt, see an example 
in Fig. I 




FIG. 4. Development of modulational instability in a 
self-focusing medium for a perturbed unstaggered nonlinear 
BW mode with Jo ~ 0.44 (a = 3, h = 0.5, and 7 = +1. 

Modulational instability of the nonlinear Bloch waves 
in a periodic medium has been earlier studied for the 
Bose-Einstein condensates in optical lattices in the 
mean-field approximation based on the Gross-Pitaevskii 
equation, which is mathematically equivalent to Eq. (|^) 
with J-"(/; x) = ^{x) + 7/. It was demonstrated that the 
unstaggered modes are always modulationally unstable in 
a self-focusing medium {x — sign7 = -|-1), and they are 
stable in a self-defocusing medium (% = — 1). It can be 
shown that the similar results are also valid for our model. 
Indeed, since the unstaggered waves are the fundamental 
modes of the self-induced periodic potential, oscillatory 
instabilities can not occur, and MI can only correspond 
to purely imaginary eigenvalues F. Such an instability 
should appear at the critical points defined by Eq. (|2J) 
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at r = 0, which are found as Q = 1, 77 + 8. Therefore, the 
range of the unstable modulation frequencies is 



7/^ 



(mm) 



< g < cos ^(77 + 3), 
< g < TT, 7] < -4, 



'4 < 7/ < -2, 



According to Eq. ([2^), for unstaggered modes we have 
XT] < — 2x, and the (in)stability results follow immedi- 
ately. We note that at small intensities (when 77 ~ —2) 
the modulational instability in a self-focusing medium 
corresponds to long- wave modulations, as illustrated in 
Fig. |. 

The staggered BW modes in a s elf- defocusing medium 
(x = —1) are also modulationally unstable p9[ |. This 
happens because the staggered waves experience effec- 
tively "normal" diffraction |^^. Such waves exist for 
77 > 2 and, similar to the case of unstaggered waves in a 
self-focusing medium, we identify the range of unstable 
frequencies corresponding to the purely imaginary eigen- 
values (ReF = 0), 

< g < cos~i(3 - 77), 2<7;<4, 
< g < TT, 77 > 4. 




FIG. 5. Modulationally unstable staggered BW modes 
(gray shading) in a self-focusing medium, shown as the inten- 
sity lo vs. the parameter a (at h = 0.5, 7 — -H). Dashed 
lines are the analytical approximations ( |2^ ) and ( p^ ) for the 
low- and high-intensity instability thresholds, respectively. 

Finally, we analyze stability of staggered BW modes in 
a self-focusing medium (x = +!)• Since such modes ex- 
ist for 77 < 2, the domain 1 < Q{T) <3-77atReF = 
does not correspond to physically possible modulation 
frequencies. However, the oscillatory instabilities (i.e. 
those with complex F) can appear due to resonances be- 
tween the modes that belong to different bands. Such 
instabilities appear in a certain region of the wave inten- 
sities in the case of shallow modulations, below a certain 
threshold value (when ah < 3.57 . . .), as shown in Fig. ||. 
We find the following asymptotic expression for the low- 
intensity instability threshold, 



while the upper boundary is given by the relation 
a ~ 47/g 'exp(7/Q '7-/4). 



(25) 



(26) 



These analytical estimates arc shown with the dashed 
lines in Fig. ^. 




FIG. 6. Region of modulational instability (gray shad- 
ing) for the model ([l9|); notations are the same as in Fig. ^. 
The result is only qualitatively similar to Fig. since the 
two-component discrete model ( p^ is asymptotically correct 
for small intensities only. 




Intensity 

FIG. 7. Unstable modulation frequencies (gray shading) 
vs. intensity Jo, for the staggered BW modes (a = 3, /i = 0.5, 
7 = -1-1). Solid line shows the parameters of the linear eigen- 
mode with the largest instability growth rate (i.e. max |Imr|). 

It is interesting to compare our results with those ob- 
tained in the fr amew ork of the continuous coupled-mode 
theory (see Sec. Ill B ), valid for the case of a narrow band 
gap (i.e. for small a and small Jp). Although the non- 
linear coupling coefficients in Eq. ( |l7| ) are different com- 
pared to a couped-mode model for shallow gratings pl| , 
the key stability result remains the same, and the oscil- 
latory instability appears above a certain critical inten- 
sity proportional to the band-gap width. In our case, 
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the band-gap width is 2a/h + 0(q;^/^), and we observe 
a good agreement with the results of the coupled-mode 
theory. However, the coupled- mode model ([l^) can not 
predict the stability region at high intensities, since in 
this region the approximation is no longer valid. 

In the limit of large a, the BW dynamics can be 
studied w ith th e help of the tight-binding approximation 
(see Sec. [II A ). Then, the effective discrete NLS equa- 
tion predicts the stability of the staggered modes 
in a self-focusing medium ||2^ . Numerical and analytical 
results confirm that our solutions arc indeed stable in the 
corresponding parameter region. 

T he tw o-component discrete model introduced in 
Sec. IIIC predicts the existence of oscillatory instabili- 
ties of the BW waves for (3 satisfying the inequality 



P2) 



id- 



^2? 



32pip2 



(/?! - P2) > 0, 



where the coefficients pi 2 are defined in Eq. (]21|). The 
corresponding instability region in the parameter space 
intensity vs. lattice depth (a) can be calculated using 
Eqs. ( p^ ) and ([20|), and the result is presented in Fig. ||. 
Since the simplified model is (asymptotically) correct 
only for small intensities (i.e. when /3 ~ there is no 
quantitative agreement between Figs. ^ and However, 
the model (|l^) does predict the key pattern of oscillatory 
MI: (i) instability appears only for a finite range of inten- 
sities when the grating depth (a) is below a critical value, 
and (ii) MI is completely suppressed for large a. Thus, 
unlike the DNLS equation, the two-component discrete 
model ( p^ predicts qualitatively all major features of MI 
in a periodic medium. 



towards the edge, 5 = tt, having there the largest in- 
stability growth rate. The corresponding modulational 
instability manifests itself through the development of 
the period-doubling modulations, as shown in Fig. H. 



V. BRIGHT SPATIAL SOLITONS 



A. Odd and Even Localized Modes 



Stationary localized modes in the form of discrete 
bright solitons can exist with the propagation constant 
inside the band gaps, when |?7| > 2. Additionally, such 
solutions can exist only if the nonlinearity and disper- 
sion sign are different, i.e. when r/x < 2. It follows from 
Eq. (^), that f3 > -{n/h)"^ and ^ > in the IR gap and 
the first BR gap (see also Fig. ||), so that the type of the 
nonlinear response is fixed by the medium characteris- 
tics, since x = sign(7). Therefore, self-focusing nonlin- 
earity can support bright solitons in the IR region (where 
77 < —2), i.e. in the conventional wave-guiding regime. In 
the case of the self-defocusing response, bright solitons 
can exist in the first BR gap, owing to the fact that the 
sign of the effective diffraction is inverted (77 > 2). In the 
latter case, the mode localization occurs in the so-called 
anti-waveguiding regime. 

Let us now consider the properties of two basic types 
of the localized modes: odd, centered at a nonlinear thin- 
film waveguide, and even, centered between the neighbor- 
ing waveguides, so that ?7|„| = x^'C^-ini-s) where s = 0, 1, 
respectively. For discrete lattices, such solutions have al- 
ready been studied in the literature (see, e.g., Ref. p3|), 
and it has been found that the mode profile is "unstag- 
gered" (i.e. C/„ > 0) if 77 < —2. On the other hand, 
Eq. (^) possesses a symmetry. 



Un 



X 



(27) 



z ^ 




FIG. 8. Development of the instability-induced pe- 
riod-doubling modulations. Initial profile corresponds to a 
slightly perturbed staggered mode with Jo — 29.87. Parame- 
ters are the same as in Fig. |^ 

We find that in a self-focusing medium the staggered 
waves {K = tt) are always stable with respect to low- 
frequency modulations, see Fig. |^. However, at larger in- 
tensities unstable frequencies appear close to the middle 
and the edge of the Brillouin zone, and they are shifted 



which means that the solutions become "staggered" at 
77 > 2. Because of this symmetry, it is sufficient to find 
localized solutions of Eq. (^) for rj < —2 and x = +1- 
The known approximate solutions give accurate results 
only in the case of highly localized modes (|ry| ^ 2) (see 
P,e3 , and references therein) and in the continuous limit 
(fryfci 2) (see, e.g., [^). In order to describe the mode 
profile for arbitrary values of 77, we introduce a new ap- 
proach based on the physical properties of localized so- 
lutions. We recall that the nonlinear localized modes 
are similar to the impurity states, which explains the 
presence of a sharp central peak (or two peaks) in the 
mode profiles at large 77. On the other hand, we have 
found that the tails of a localized mode are always quite 
smooth, both in the continuous limit and highly-discrete 
case. Based on these facts, we construct the approximate 
solutions by matching the mode tails with the central im- 
purity node(s). 
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First, we have to find an approximation for the mode 
tails. Since the tail is smooth, its profile can be well ap- 
proximated by the continuous equations. The simplest 
yet effective approach is to choose the model coefficients 
to match the discrete solutions at the beginning of the 
tail, which we define as the zero concavity point, and 
in the linear limit corresponding to the far-field asymp- 
totics. Then, after simple calculations, we obtain an ap- 
proximate continuous equation for the mode tails, 



-XU- 

where A = — (r/ + 2) 
solution of Eq. (gSj) 
totics has the form 



A d^U 
drfi 



= 0, 



(28) 



> 0, and p = cosh"^(l + A/2). The 
with the vanishing far-field asymp- 



(29) 



C/(n; Us) = v2Asech(p(n + Us)), 



where the free parameter defines the shift. Note that 
in the limit A — > we have p VA, and the results of the 
conventional continuous approximation ||l| are recovered. 
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1 



1 2 3 4 5 

FIG. 9. Dependence of (a) the shift parameter and 
(b) peak normalized amplitude [/o in Eqs. (^o|), ( ^l|) on the 
parameter A for odd and even bright localized modes. Dashed 
lines in (b) — numerically calculated values. 

Second, we have to construct the full approximate so- 
lution as a combination of two tails. Such a state is 
supported by the central node in odd modes (J7o), and 
by two nodes in the case of even topology {U-\ = xC/o)- 
Since the profiles are symmetric, we have to calculate 
the field structure only for n > 0, and the discrete tail 
profiles for n > 1 can be approximated as: 



(30) 



where the function U{n; Ug) is given by Eq. (|29|). In 
order to determine the unknown parameters, the peak 
amplitude Uq and shift Ug, we should solve the original 



discrete equations (||) at the soliton peak {n = 0) and at 
the neighboring site (n = 1): 



-(2 - s + \)Uo + (2 - s)Ui + ^ 0, 
^(2 + A)[/i + J7o + t/2 + t/i -0, 



(31) 



where we took into account the symmetry properties of 
odd and even modes. We find that for all A > there 
exists a solution that belongs to a smooth branch, start- 
ing with Hg = s/2 in the continuous limit (at A — > 0), see 
Fig. |9|(a) . We compare the analytical approximation with 
the exact numerical solution of the original model (|^) 
for the peak amplitude Uq, and find that an error does 
not exceed 1.5% for odd and 0.8% for even modes. As 
a matter of fact, the corresponding curves are indistin- 
guishable in Fig. ||(b). The actual mode profiles are 
also adequately represented, see examples in Figs. p^(a,b) 
and |ll|(a,b). Therefore, the suggested analytical proce- 
dure allows us to obtain extremely accurate approximate 
analytical solutions in the whole parameter range, in- 
cluding the extreme cases of the continuous (A — s- 0) and 
anti-continuous (A +00) limits. 






x+h/2 

FIG. 10. Top: power vs. propagation constant for 
odd (black) and even (dashed gray) localized modes in a 
self-focusing (7 — +1) regime: solid — stable, dashed — 
unstable. Gray shading marks the transmission band. Bot- 
tom: profiles of the localized modes corresponding to the 
marked points (a) and (b) in the top plot; black dots — an- 
alytical approximation for the node amplitudes. The lattice 
parameters are the same as in Fig. ^ 

Our linear stability analysis reveals that even modes 
are always unstable with respect to a translational shift 
along the x axis. On the other hand, odd modes are 
always stable in the self-focusing regime (see Fig. |l^, 
top), but can exhibit oscillatory instabilities in the self- 
defocusing case when the power exceeds a certain critical 
value (see Fig. |ll|, top). At this point, an eigenmode of 
the linearized problem resonates with the band-gap edge, 
the value (/3 — Re F) moves inside the band, and non-zero 
imaginary part of the eigenvalue appears, as illustrated 
by an example in Fig. 02. Such an instability scenario is 
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similar to one earlier identified for gap solitons |2^, and 
also for the modes localized at a single nonlinear layer 
in a linear periodic structure [|o|. The latter example 
demonstrates a deep similarity between the periodic sys- 
tems with localized and distributed nonlinearities. 
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FIG. 11. Top: power vs. propagation constant in the 
self-defocusing (7 = —1) regime. Dotted — oscillatory unsta- 
ble modes; other notations are the same as in Fig. hol 
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FIG. 12. Example of a resonance between an eigenmode 
of the linear eigenvalue problem and a bang-gap edge of the 
continuous spectrum that leads to an oscillatory instability of 
the odd localized mode (parameters correspond to Fig. hll). 



Development of both types of instability is demon- 
strated in Figs. |l^(a,b). Transformation of an even mode 
into an odd counterpart due to a symmetry-breaking in- 
stability is shown in Fig. ^(a). The effect of oscillatory 
instability on an odd mode is quite different: strong ra- 
diation is emitted due to a resonant coupling with linear 
waves outside the band gap, see Fig. |l|(b). 




FIG. 13. Instability development for the Bragg-type 
localized modes in a self-defocusing medium: (a) symme- 
try-breaking instability of an even mode {P ~ 3.23); (b) os- 
cillatory instability of an odd mode (P ~ 23.4). Parameters 
correspond to Fig. |l^. 



B. Soliton Bound States 



"Twisted" Modes 



Due to a periodic modulation of the medium refractive 
index, solitons can form bound states [P6| . In particular, 
the so-called "twisted" localized mode p7| i s a combina- 
tion of two out-of-phase bright solitons p6| . Such solu- 
tions do not have their continuous counterparts, and they 
can only exist when the discreteness effects are strong, i.e. 
for |?7| > rjcT- Properties of the twisted modes depend on 
the separation between the modes forming a bound state. 
We consider the cases of two lowest-order solutions of 
(i) "even" type with zero nodes (m = 0) in-between the 
peaks, and (ii) "odd" type with one node (m = 1) at 
the middle, with Uq = 0. The corresponding symmetry 
properties are C/|„|_|_,„ = — X™^^t^-|ri|-i- Additionally, 
Eq. (^) also holds, so that we only have to construct 
solutions for x = +1- Then, the and the soliton tails at 
n > m are approximated as: 



Un = U{n ~ m; Ug), 



(32) 



where U{n; Ug) is given by Eq. (p9|). The matching con- 
ditions are: 
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-(3 - m + A)C/„ + Urn+i + Ui = 0, 

-(2 + X)Ura+l + U,n + U.m+2 + 

where, as before, A = — (t; + 2) > 0. 
0.6 



0, 



(33) 




FIG. 14. Dependence of (a) the shift parameter Us and 
(b) the peak normalized amphtude Um in Eqs. ( p^ ) and ( p^ ) 
on the parameter A, for odd and even twisted locahzed modes. 
Dashed hnes in (b) — numerically calculated values. 



the best of our knowledge, none of the previously devel- 
oped analytical approximations could predict the regions 
of existence for highly discrete twisted modes (with small 
to). Moreover, the approximate solution describes very 
accurately the profiles of the twisted modes, see exam- 
ples in Figs. |lj(a,b) and Figs. |l^(a,b). In particular, 
a relative error for the peak amplitude Um is less than 
0.5%, so that in Fig. ^(b) the analytical and numerical 
dependencies practically coinside. 
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FIG. 16. Top: power vs. propagation constant for the odd 
twisted localized waves in a self-defocusing (7 = — 1) regime. 
Notations are the same as in Figs. |l^ and|ll} 






x+h/2 

FIG. 15. Top: power vs. propagation constant for 
odd (black) and even (gray) twisted localized waves in a 
self-focusing (7 = +1) regime. Notations are the same as 
in Figs. |l| and 0. 

We determine solution of Eqs. (|33| ) starting with the 
highly localized modes earlier described anti-continuous 
limit (A 3> 1) [§7|, and then gradually decrease the pa- 
rameter A. We find that solution exists for A > Ad- > 0, 
and disappears when rig (Ad ) = 0. Substituting this con- 
dition into Eq. (|33|), we determine the approximate crit- 
ical parameter values, rj^im = 0) ~ 3.32 and rjcrim = 
1) ~ 2.95. The analytically determined existence regions 
agree very well with numerical results (within 1%). To 



In the self-focusing regime, stability properties of the 
twisted modes in the IR gap can be similar to those ear- 
lier identified in the framework of a DNLS model ||2^. In 
the example shown in Fig. 15, the modes are stable at 



larger values of the propagation constant, and they be- 
come oscillatory unstable closer to the boundary of the 
existence region. Quite importantly, the stability region 
is much wider in the case of odd twisted modes due to a 
larger separation between the individual solitons of the 
bound state. 



The characteristics of the twisted modes in the BR gap 
can differ substantially from the previous case. First, 
the value of r/ is limited from above (2 < 77 < ?7max) 
and, therefore, some families of the twisted modes with 
TO < TOcr may not exist. For example, for the medium 
parameters corresponding to Fig. |^(a), we have rjcr{m — 
1) < ('7max — 3.18) < 77cr(™ — 0), SO that TOcr = 1. Un- 
der these conditions, even modes with m = cannot ex- 
ist in the BR regime. On the other hand, the odd modes 
with m = 1 can exist, and they are stable in a wide pa- 
rameter region, see Fig. |l6|. Development of oscillatory 
instabilities for IR and BR twisted modes is illustrated 
in Figs. |l^(a) and |l^(b), respectively. In both the cases, 
we observe an exponential increase of periodic amplitude 
modulations and the emission of radiation waves. 
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(a) 





FIG. 17. Instability scenarios for twisted localized modes: 
(a) an even mode in a self-focusing medium (P ~ 1.22): (b) an 
odd mode in a self-defocusing regime (P ~ 4.45). Parameters 
correspond to Figs, and hfJ, respectively. 



VI. DARK SPATIAL SOLITONS 

Similar to the continuous NLS equation with self- 
defocusing nonhnearity |^ or the DNLS equation [p9| , 
our model can support dark solitons — localized modes 
on the Bloch-wave background. However, dark station- 
ary localized modes in a periodic medium can exist for 
both signs of nonlinearity. To be specific, let us consider 
the case of a background corresponding to the Bloch- 
wave solutions with K — 0,Tr introduced in Sec. IV. 
Then, dark-mode solutions can appear at the band-gap 
edge where = 2x = 2sign7, since in such a case nonlin- 
ear and dispersion terms have the same signs 

Similar to the case of bright solitons discussed above, 
two basic types of dark spatial solitons can be identi- 
fied, namely, odd localized modes centered at a nonlin- 
ear thin-film waveguide, and even localized modes cen- 
tered between the neighboring thin-film waveguides. All 



such modes satisfy the symmetry condition, Ui 



\n\+s 



~{~x)^ U^\n\-i, where s = 0, 1 for even and odd 
modes, respectively. The BW background is unstaggered 
if X = —1, and it is staggered for x = -1-1; the correspond- 
ing solutions can be constructed with the help of a sym- 



metry transformation, Un — > (— l)"C/„. However, the sta- 
bility properties of these two types of localized states can 
be quit e different. Indeed, it has been demonstrated in 
Sec. f\1that in a self- focusing medium (x = +1) the stag- 
gered background can become unstable. On the contrary, 
the unstaggered background is always stable if x = — 1- 

In order to find the approximate analytical solutions, 
we consider the case x = ^1, with no lack of generality 
[since solutions with x = +1 can be obtained by ap- 
plying the symmetry transformation (p^)]. In this case, 
the far-field asymptotics for solutions of Eq. (||) close 
to the background level can be found as {Uoo — Un) — 
eP", where Uca = \/A is the background amplitude, 
p ~ cosh~^(l + A) is the localization parameter, and 
A = 7^ -f 2 > 0. Then, we obtain an approximate contin- 
uous equation for the nonlinear mode tails by matching 
the asymptotic solution at large n, 



XU 



2Xd^U 



-u-" = 0. 



(34) 



The corresponding dark-soliton solution has the form 



U{n; Us) = VA tanh(p(n -I- ns)/2). 



(35) 



Note that in the limit A ^ we have p ^ v 2A, and the 
results of the conventional continuous approximation are 
recovered. 

Similar to the discrete bright solitons, a localized so- 
lution can be constructed by matching the soliton tails 
defined by Eq. ( pO| ) . The corresponding matching condi- 
tions, 



(A - 2 - s)Uo + (1 - s)Ui - = 0, 
(A - 2)Ui + U0 + U2- Uf = 0, 



(36) 



are used to determine the shift parameter and ampli- 
tude Uq. We have Uq = for odd modes, due to their 
symmetry properties. Dependencies of the shift param- 
eter Us on A are presented in Fig. 18(a). The soliton 
amplitudes at the central sites are shown in Fig. p^(b), 
where we observe again an excellent agreement between 
the approximate analytical and numerical solutions (the 
corresponding errors do not exceed 2% for odd and 1% 
for even modes). 
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FIG. 18. Dependence of (a) the shift parameter Us and 
(b) the normalized amphtudes Us in Eqs. (|3o|),(|3^), and (^) 
on the parameter A for odd (s — 1) and even (s = 0) twisted 
locahzed modes. Dashed hnes in (b) — numerically calculated 
values. 

Two types of dark spatial solitons in our model are 
shown for both staggered and unstaggered BW back- 
grounds in Figs. |l^ and respectively. We characterize 
the family of dark solitons by the complimentary power 
defined as 

/-\-nh 
{\u{x + 2n/i)p - |w(a;)p) dx, 
-nh 



where n is integer. The localized solutions shown in 
Fig. ^ are similar to those found earlier in Ref. |^ in 
the context of the superflow dynamics on a periodic po- 
tential. 




FIG. 20. Top: complementary power vs. propagation con- 
stant in the self-defocusing (7 — —1) regime. Notations are 
the same as in Fig. [13. 
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FIG. 19. Top: complementary power vs. propagation con- 
stant for odd (black) and even (gray) dark localized solitons 
in a self-focusing (7 = -1-1) regime. Notations are the same as 
in Fig. hol but (in)stability regions are not indicated. 
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FIG. 21. Propagation dynamics of (a) even and (b) odd 
dark localized modes in a self- focusing medium. Initial pro- 
files correspond to slightly perturbed stationary solutions at 
13 = 20, other parameters correspond to Fig. 



Numerical study of the propagation dynamics demon- 
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strates that, similar to the case of bright soHtons, even 
dark-soliton modes are unstable with respect to asym- 
metric perturbations, see an example in Fig. ^(a). On 
the other hand, odd modes can propagate in a stable (or 
weakly unstable) manner, as illustrated in Fig. |l|(b). We 
note however that dark solitons can exhibit oscillatory 
instabilities close to the continuum limit (at small 
intensities), but a detailed analysis of the dark- mode sta- 
bility is beyond the scope of the present paper. 



VII. CONCLUSION 

In the framework of a simplified model of a nonlin- 
ear layered medium that describes the so-called Dirac- 
comb nonlinear waveguide array, we have analyzed spa- 
tial optical solitons in the form of bright, dark, and 
"twisted" localized modes. In general, such solitons are 
of two types, i.e. they are either (i) nonlinear guided 
waves localized due to the total internal reflection or 
(ii) the Bragg-type localized modes existing in the for- 
bidden transmission gaps, gap solitons. We have ana- 
lyzed the existence and stability of the nonlinear localized 
modes of both types and described also modulational in- 
stability of extended modes induced by a periodic change 
of the medium refractive index. Additionally, we have 
discussed both similarities and differences with the mod- 
els described by the DNLS equation, derived in the fre- 
quently used tight-binding approximation, and with the 
results of the coupled-mode theory, which are valid for 
a shallow modulation and a narrow gap in the transmis- 
sion spectrum. We believe our analysis and results may 
be useful for other fields, such as the nonlinear dynamics 
of the Bose-Einstein condensates in optical lattices (see, 
e.g., Ref. i). 
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